clear
capture log close#delimit ;set more off;version 11.0;

*************************************************************;
*	Replication files for Conrad & Ritter (2013)	    *;
*************************************************************;
use "ConradRitterAnnualData19June12.dta", clear;
sort cowcode year;
**Generating Interactions**;
gen JEJS=linzerstatonJI*jobsecurity_cheibub;

capture program drop biselect;
program define biselect;
    args lnf theta1 theta2 theta3 zrhoS zrhoN;

	tempvar rhoN rhoS;
	qui gen double `rhoN' = (exp(2*`zrhoN')-1) / (exp(2*`zrhoN')+1); 
	qui gen double `rhoS' = (exp(2*`zrhoS')-1) / (exp(2*`zrhoS')+1);

quietly replace `lnf' = ln(($ML_y1==1)*binorm(`theta1',$ML_y2*`theta2', $ML_y2*`rhoS')+($ML_y1==0)*binorm(-`theta1',$ML_y3*`theta3', -$ML_y3*`rhoN'));
          
end;

set matsize 150;
capture drop y1;
capture drop y2;
capture drop y3;
capture drop y1n;

gen y1=catr;
gen y1n=-y1;
gen y2=2*CIRIhightort -1;
gen y3=y2;
	
local ivar1 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS IO;
local ivar2 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;
local ivar3 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;

prob y1 `ivar1';
local ll1=e(ll);
qui matrix i1 = e(b);
qui matrix coleq i1 = eq1:;

prob CIRIhightort `ivar2'; 
local ll2=e(ll);
qui matrix i2 = e(b);
qui matrix coleq i2 = eq2:;

prob CIRIhightort `ivar3' if y1==0;
local ll3=e(ll);
qui matrix i3 = e(b);
qui matrix coleq i3 = eq3:;
 
heckprob CIRIhightort `ivar2', select (y1 = `ivar1'); 
qui matrix i4 = e(b);
qui matrix rS=i4[1,"athrho:_cons"];
qui matrix coleq rS= eq4:;

heckprob CIRIhightort `ivar3', select (y1n = `ivar1');
qui matrix i5 = e(b);
qui matrix rN=-i5[1,"athrho:_cons"];
qui matrix coleq rN= eq5:;

qui matrix initvec = i1,i2,i3,rS,rN;

ml model lf biselect (y1 = `ivar1') (y2 = `ivar2') (y3 = `ivar3') () ();
ml search;
ml maximize, difficult;

est store est1;

disp (exp(2*[eq4]_cons)-1) / (exp(2*[eq4]_cons)+1);
disp (exp(2*[eq5]_cons)-1) / (exp(2*[eq5]_cons)+1);

*************************************************************;

capture program drop biselect;
program define biselect;
version 7;
         args lnf theta1 theta2 theta3  ;

quietly replace `lnf' = ln(($ML_y1==1)*binorm(`theta1',$ML_y2*`theta2', $ML_y2*0)+($ML_y1==0)*binorm(-`theta1',$ML_y3*`theta3', -$ML_y3*0));
end;

set matsize 150;
capture drop y1;
capture drop y2;
capture drop y3;
capture drop y1n;

gen y1=catr;
gen y1n=-y1;
gen y2=2*CIRIhightort-1;
gen y3=y2;

local ivar1 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS IO;
local ivar2 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;
local ivar3 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;

prob y1 `ivar1' ;
local ll1=e(ll);
qui matrix i1 = e(b);
qui matrix coleq i1 = eq1:;

prob CIRIhightort `ivar2' if y1==1;
local ll2=e(ll);
qui matrix i2 = e(b);
qui matrix coleq i2 = eq2:;
 
prob CIRIhightort `ivar3' if y1==0;
local ll3=e(ll);
qui matrix i3 = e(b);
qui matrix coleq i3 = eq3:;

heckprob CIRIhightort `ivar2', select (y1 = `ivar1') ;
qui matrix i4 = e(b);
qui matrix rS=i4[1,"athrho:_cons"];
qui matrix coleq rS= eq4:;

heckprob CIRIhightort `ivar3', select (y1n = `ivar1');
qui matrix i5 = e(b);
qui matrix rN=-i5[1,"athrho:_cons"];
qui matrix coleq rN= eq5:;

qui matrix initvec = i1,i2,i3;

ml model lf biselect (y1 = `ivar1') (y2 = `ivar2') (y3 = `ivar3');
ml init initvec, skip;
ml search;
ml maximize, difficult;

*CALCULATE THE PROBABILITY THAT JOINT EFFECT OF RHO-S AND RHO-N = 0;
est store est2;
lrtest est1 est2;

*************************************************************;

capture program drop biselect;
program define biselect;
version 7;
         args lnf theta1 theta2 theta3 zrhoN ;

	tempvar rhoN rhoS ;
	qui gen double `rhoN' = (exp(2*`zrhoN')-1) / (exp(2*`zrhoN')+1) ;

      quietly replace `lnf' = ln(($ML_y1==1)*binorm(`theta1',$ML_y2*`theta2', $ML_y2*0)+($ML_y1==0)*binorm(-`theta1',$ML_y3*`theta3', -$ML_y3*`rhoN'));
end;

set matsize 150;
capture drop y1;
capture drop y2;
capture drop y3;
capture drop y1n;

gen y1=catr;
gen y1n=-y1;
gen y2=2*CIRIhightort -1;
gen y3=y2;

local ivar1 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS IO;
local ivar2 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;
local ivar3 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;

prob y1 `ivar1' ;
local ll1=e(ll);
qui matrix i1 = e(b);
qui matrix coleq i1 = eq1:;

prob CIRIhightort `ivar2' if y1==1;
local ll2=e(ll);
qui matrix i2 = e(b);
qui matrix coleq i2 = eq2:;
 
prob CIRIhightort `ivar3' if y1==0;
local ll3=e(ll);
qui matrix i3 = e(b);
qui matrix coleq i3 = eq3:;

heckprob CIRIhightort `ivar2', select (y1 = `ivar1') ;
qui matrix i4 = e(b);
qui matrix rS=i4[1,"athrho:_cons"];
qui matrix coleq rS= eq4:;

heckprob CIRIhightort `ivar3', select (y1n = `ivar1');
qui matrix i5 = e(b);
qui matrix rN=-i5[1,"athrho:_cons"];
qui matrix coleq rN= eq5:;

qui matrix initvec = i1,i2,i3,rN;

ml model lf biselect (y1 = `ivar1') (y2 = `ivar2') (y3 = `ivar3') ();
ml init initvec, skip;
ml search;
ml maximize, difficult;

*OBTAIN P-VALUE ASSOCIATED WITH RHO-S:;
est store est3;
lrtest est1 est3;

*************************************************************;

capture program drop biselect;
program define biselect;
version 7;
         args lnf theta1 theta2 theta3 zrhoS ;

	tempvar rhoN rhoS ;
	qui gen double `rhoS' = (exp(2*`zrhoS')-1) / (exp(2*`zrhoS')+1);

      quietly replace `lnf' = ln(($ML_y1==1)*binorm(`theta1',$ML_y2*`theta2', $ML_y2*`rhoS')+($ML_y1==0)*binorm(-`theta1',$ML_y3*`theta3', -$ML_y3*0));
end;

set matsize 150;
capture drop y1;
capture drop y2;
capture drop y3;
capture drop y1n;

gen y1=catr;
gen y1n=-y1;
gen y2=2*CIRIhightort -1;
gen y3=y2;

local ivar1 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS IO;
local ivar2 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;
local ivar3 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;

prob y1 `ivar1';
local ll1=e(ll);
qui matrix i1 = e(b);
qui matrix coleq i1 = eq1:;

prob CIRIhightort `ivar2' if y1==1;
local ll2=e(ll);
qui matrix i2 = e(b);
qui matrix coleq i2 = eq2:;
 
prob CIRIhightort `ivar3' if y1==0;
local ll3=e(ll);
qui matrix i3 = e(b);
qui matrix coleq i3 = eq3:;

heckprob CIRIhightort `ivar2', select (y1 = `ivar1') ;
qui matrix i4 = e(b);
qui matrix rS=i4[1,"athrho:_cons"];
qui matrix coleq rS= eq4:;

heckprob CIRIhightort `ivar3', select (y1n = `ivar1');
qui matrix i5 = e(b);
qui matrix rN=-i5[1,"athrho:_cons"];
qui matrix coleq rN= eq5:;

qui matrix initvec = i1,i2,i3,rS;

ml model lf biselect (y1 = `ivar1') (y2 = `ivar2') (y3 = `ivar3') ();
ml init initvec, skip;
ml search;
ml maximize, difficult;

*OBTAIN P-VALUE ASSOCIATED WITH RHO-N;
est store est4;
lrtest est1 est4;
exit;